Introduction à la modélisation statistique bayésienne

Planning

Cours n°01 : Introduction à l’inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux
Cours n°09 : Modèles multi-niveaux généralisés
Cours n°10 : Data Hackathon

\[\newcommand\given[1][]{\:#1\vert\:}\]

Le problème avec la distribution postérieure

\[ \color{purple}{p(\mu, \sigma \given h)} = \frac{\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \given \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i} \given \mu, \sigma)\mathrm{Normal}(\mu \given 178, 20)\mathrm{Uniform}(\sigma \given 0, 50) \mathrm{d} \mu \mathrm{d} \sigma}} \]

Petit problème : \(p(\text{data})\) s’obtient en calculant la somme (pour des variables discrètes) ou l’intégrale (pour des variables continues) de la densité conjointe \(p(\text{data}, \theta)\) sur toutes les valeurs possibles de \(\theta\). Cela se complique lorsque le modèle comprend plusieurs paramètres.

Le problème avec la distribution postérieure

Rappels Cours n°02

Trois méthodes pour résoudre (contourner) ce problème :

• La distribution a priori est un prior conjugué de la fonction de vraisemblance (e.g., modèle Beta-Binomial). Dans ce cas, il existe une solution analytique (i.e., qu’on peut calculer de manière exacte) pour la distribution postérieure.

• Autrement, pour des modèles simples, on peut utiliser la méthode par grille. On calcule la valeur exacte de la probabilité postérieure en un nombre fini de points dans l’espace des paramètres.

• Pour les modèles plus complexes, explorer tout l’espace des paramètres n’est pas tractable. On va plutôt échantillonner intelligemment un grand nombre de points dans l’espace des paramètres.

Objectifs du cours

\(\longrightarrow~\) Présenter le principe de base de l’échantillonnage : Markov Chain Monte Carlo

\(\longrightarrow~\) Présenter deux algorithmes (Metropolis-Hastings et HMC)

\(\longrightarrow~\) Montrer les forces mais aussi les faiblesses de ces méthodes

\(\longrightarrow~\) Donner des outils de contrôle sur ces méthodes

\(\longrightarrow~\) Appliquer ces méthodes à un cas simple

Markov Chain Monte Carlo

• Markov chain Monte Carlo
  \(\longrightarrow~\) Échantillonnage aléatoire
  \(\longrightarrow~\) Le résultat est un ensemble de valeurs du paramètre

• Markov chain Monte Carlo
  \(\longrightarrow~\) Les valeurs sont générées sous forme de séquences (liaison de dépendance)
  \(\longrightarrow~\) Indice temporel pour identifier la place dans la chaîne
  \(\longrightarrow~\) Le résultat est de la forme : \(\theta^1, \theta^2, \theta^3, \dots, \theta^t\)

• Markov chain Monte Carlo
  \(\longrightarrow~\) La valeur de paramètre générée ne dépend que de la valeur du paramètre précédent \(\Pr(\theta^{t+1} \given \theta^{t}, \theta^{t-1}, \ldots, \theta^{1}) = \Pr(\theta^{t + 1} \given \theta^{t})\)

Méthodes Monte Carlo

Le terme de méthode de Monte-Carlo désigne une famille d’algorithmes visant à calculer (ou approcher) une valeur numérique en utilisant des procédés aléatoires, c’est-à-dire des techniques probabilistes. Cette méthode a été formalisée en 1947 par Nicholas Metropolis, et publiée pour la première fois en 1949 dans un article co-écrit avec Stanislaw Ulam.

Méthodes Monte Carlo : Estimation de \(\pi\)

Soit un point \(M\) de coordonnées \((x, y)\), où \(0 < x < 1\) et \(0 < y < 1\). On tire aléatoirement les valeurs de \(x\) et \(y\) entre \(0\) et \(1\) suivant une loi uniforme. Le point \(M\) appartient au disque de centre \((0, 0)\) de rayon \(R = 1\) si et seulement si \(x^{2} + y^{2} \leqslant 1\). La probabilité que le point \(M\) appartienne au disque est \(\pi / 4\) puisque le quart de disque est de surface \(\sigma = \frac{\pi R^{2}}{4} = \frac{\pi}{4}\), et le carré qui le contient est de surface \(S = R^{2} = 1\). Si la loi de probabilité du tirage de point est uniforme, la probabilité de tomber dans le quart de disque vaut \(\frac{\sigma}{S} = \frac{\pi}{4}\). En faisant le rapport du nombre de points dans le disque au nombre de tirages, on obtient une approximation du nombre \(\pi / 4\) si le nombre de tirages est grand.

Méthodes Monte Carlo : Estimation de \(\pi\)

trials <- 1e5
radius <- 1
x <- runif(n = trials, min = 0, max = radius)
y <- runif(n = trials, min = 0, max = radius)
distance <- sqrt(x^2 + y^2)
inside <- distance < radius
pi_estimate <- 4 * sum(inside) / trials

Méthodes Monte Carlo

Autre exemple : déterminer la superficie d’un lac ou encore déterminer le maximum d’une fonction (optimisation) via simulated annealing ou recuit simulé (Wikipedia).

Méthodes Monte Carlo

Monte Carlo désigne une famille d’algorithmes qui ont pour but d’approcher des valeurs numériques à partir de procédés aléatoires. Pourrait-on s’en servir pour obtenir une approximation de la distribution postérieure ?

Influence de la constante de normalisation

Méthodes Monte Carlo : Exemple

Considérons un exemple simple : Soit un paramètre \(\theta\) avec 7 valeurs possibles et la fonction de répartition suivante.

Méthodes Monte Carlo : Exemple

Estimation de cette distribution par tirage aléatoire : Cela revient à tirer aléatoirement un grand nombre de points “au hasard” parmi ces 28 cases (comme pour le calcul de \(\pi\)) !

Méthodes Monte Carlo : Exemple

niter <- 100
theta <- 1:7
ptheta <- theta
samples <- sample(x = theta, prob = ptheta, size = niter, replace = TRUE)

• La distribution des échantillons obtenus converge vers la “vraie” distribution
• Mais, nécessite généralement beaucoup d’échantillons…
• Aucun contrôle sur la vitesse de convergence…
• Et si on abandonnait l’échantillonnage indépendant ?

Algorithme Metropolis

Cet algorithme a été présenté pour la première fois en 1953 par Nicholas Metropolis, Arianna W. Rosenbluth, Marshall Rosenbluth, Augusta H. Teller, et Edward Teller. Le problème des algorithmes Monte-Carlo n’est pas la convergence, mais la vitesse à laquelle la méthode converge. Pour augmenter la vitesse de convergence, il faudrait faciliter l’accès aux valeurs de paramètres les plus représentées.

Algorithme Metropolis

Sélectionner un point de départ (on peut sélectionner n’importe quelle valeur).

Algorithme Metropolis

Faire une proposition de déplacement centrée sur la valeur courante de \(\theta\).

Algorithme Metropolis

Calculer la probabilité d’accepter le déplacement :

\[\Pr_{\text{move}} = \text{min} \left(\frac{\Pr(\theta_{\text{proposed}})}{\Pr(\theta_{\text{current}})}, 1 \right)\]

Algorithme Metropolis

La position calculée devient la nouvelle position de départ et on répète l’algorithme.

Algorithme Metropolis

metropolis <- function (niter = 1e2, startval = 3) {
    
    x <- rep(0, niter) # initialises the chain vector
    x[1] <- startval # defines the starting value
    
    for (i in 2:niter) {
        
        current <- x[i - 1] # current value of the parameter
        proposal <- current + sample(c(-1, 1), size = 1)
        # making sure proposals are within the possible range
        if (proposal < 1) proposal <- 1
        if (proposal > 7) proposal <- 7
        prob_move <- min(1, proposal / current)
        x[i] <- sample(c(current, proposal), size = 1, prob = c(1 - prob_move, prob_move) )
        
    }
    
    return (x)
    
}

Méthodes Monte Carlo vs. Algorithme Metropolis

Algorithme Metropolis

Application au lancer de pièce (cas continu)

\(~\bullet~\) La fonction de vraisemblance est donnée par : \(\color{orangered}{p(y \given \theta, n) = \theta^y(1 - \theta)^{(n - y)}}\)
\(~\bullet~\) Le prior est donné par : \(\color{steelblue}{p(\theta \given a, b) \propto \theta^{(a - 1)}(1 - \theta)^{(b - 1)}}\)
\(~\bullet~\) Le paramètre que l’on cherche à estimer prend ses valeurs dans l’intervalle \(\left[0, 1 \right]\)

Algorithme Metropolis

Problème n°2 : Quelle probabilité utiliser pour accepter ou refuser le déplacement ? Nous utilisons le produit de la vraisemblance et du prior : \(\color{orangered}{\theta^y(1 - \theta)^{(n - y)}}\color{steelblue}{\theta^{(a - 1)}(1 - \theta)^{(b - 1)}}\)

Algorithme Metropolis

\(\longrightarrow~\) Sélectionner un point de départ
\(~\bullet~\) Il faut choisir \(\theta \in \left[0,1\right]\)
\(~\bullet~\) Seule contrainte \(\Pr(\theta_{initial}) \ne 0\)

Algorithme Metropolis

Algorithme Metropolis

Le choix de sigma dans la proposition de déplacement

Deux indices permettent d’évaluer la qualité de l’échantillonnage :

\(\rightarrow\) Le rapport entre le nombre de déplacements proposés et le nombre de déplacements acceptés

\(\rightarrow\) L’effective sample size (i.e., le nombre de déplacements qui ne sont pas corrélés avec les précédents)

Algorithme Metropolis

Algorithme Metropolis

Le choix de sigma dans la proposition de déplacement

\(\rightarrow~\) Toutes les propositions de déplacement (ou presque) sont acceptées

\(\rightarrow~\) Peu de valeurs effectives

Il faut beaucoup d’itérations pour avoir un résultat satisfaisant…

Algorithme Metropolis

Algorithme Metropolis

Le choix de sigma dans la proposition de déplacement

\(\rightarrow\) Les propositions de déplacement sont rarement acceptées

\(\rightarrow\) Peu de valeurs effectives…

Il faut beaucoup d’itérations pour obtenir un résultat satisfaisant…

Algorithme Metropolis-Hastings

metropolis_hastings <- function (niter = 1e2, startval = 0.5) {
    
    x <- rep(0, niter) # initialises the chain vector
    x[1] <- startval # defines the starting value
    
    for (i in 2:niter) {
        
        current <- x[i - 1] # current value of the parameter
        current_plaus <- dbeta(current, 2, 3) * dbinom(1, 2, current)
        # proposal <- runif(n = 1, min = 0, max = 1) # proposed value of the parameter
        proposal <- rnorm(n = 1, mean = current, sd = 0.1) # proposed value of the parameter
        # making sure proposals are within the possible range
        if (proposal < 0) proposal <- 0
        if (proposal > 1) proposal <- 1
        proposal_plaus <- dbeta(proposal, 2, 3) * dbinom(1, 2, proposal)
        alpha <- min(1, proposal_plaus / current_plaus) # moving probability ratio
        x[i] <- sample(c(current, proposal), size = 1, prob = c(1 - alpha, alpha) )
        
    }
    
    return (x)
    
}

Algorithme Metropolis-Hastings

z1 <- metropolis_hastings(niter = 1e4, startval = 0.5)
z2 <- metropolis_hastings(niter = 1e4, startval = 0.5)

data.frame(z1 = z1, z2 = z2) %>%
  mutate(sample = 1:1e4) %>%
  pivot_longer(cols = z1:z2) %>%
  ggplot(aes(x = sample, y = value, colour = name) ) +
  geom_line(show.legend = FALSE) +
  labs(x = "Nombre d'itérations", y = expression(theta) )

Algorithme Metropolis-Hastings

data.frame(z1 = z1, z2 = z2) %>%
  pivot_longer(cols = z1:z2) %>%
  rownames_to_column() %>%
  mutate(rowname = as.numeric(rowname) ) %>%
  ggplot(aes(x = value) ) +
  geom_histogram(aes(y = ..density..), color = "white", alpha = 0.8) +
  stat_function(fun = dbeta, args = list(3, 4), color = posterior_color, size = 1) +
  facet_wrap(~name) +
  labs(x = expression(theta), y = "Nombre d'échantillons")

Algorithme Metropolis-Hastings

Algorithme Hamiltonian Monte Carlo

Les algorithmes Metropolis et Metropolis-Hastings (ou Gibbs) ont de mauvaises performances lorsque les paramètres du modèle sont fortement corrêlés. L’algorithme Hamiltonian Monte Carlo repose sur l’algorithme Metropolis mais résout ces problème en utilisant la géométrie de l’espace postérieur. On utilise l’opérateur hamiltonien (hamiltonians) qui représente l’énergie totale d’un système. Cette énergie se décompose en l’énergie potentielle (qui dépend de la position dans l’espace des paramètres \(\theta\)) et son énergie cinétique, qui dépend de son momentum (\(m\)) :

\[ H(\theta, m) = \underbrace{U(\theta)}_{\text{énergie potentielle}} + \underbrace{KE(m)}_{\text{énergie cinétique}} \]

Algorithme Hamiltonian Monte Carlo

Si on “suit” simplement la gravité de notre position dans l’espace postérieur (imaginons que nous suivions la trajectoire d’une luge ou d’une bille dans cet espace), on arrivera à un minimum de cet espace, mais on veut explorer tout l’espace proportionnellement à la densité postérieure, pas seulement les endroits de forte densité…

Algorithme Hamiltonian Monte Carlo

Algorithme Hamiltonian Monte Carlo

\(\longrightarrow~\) Sélectionner un point de départ \(\theta_{0}\) : On peut sélectionner n’importe quelle valeur de \(\theta\) dans l’espace des paramètres

Influence de la durée de déplacement…

Influence de la variabilité du moment initial…

Algorithme Hamiltonian Monte Carlo

Algorithme Hamiltonian Monte Carlo

• L’algorithme HMC est particulièrement utile lorsque le modèle contient un grand nombre de paramètres et quand les paramètres sont corrélés entre eux.

• Il nécessite un paramétrage assez fin. Des logiciels proposent des méthodes clés en main très efficaces.
  

• Le package brms utilise une variation de cet algorithme (l’algorithme NUTS) qui permet d’éviter les trajectoires cycliques.

Évaluation des MCMC

Ces méthodes peuvent ne pas converger vers la “vraie” distribution postérieure, en raison du temps de calcul limité, du paramétrage de certains hyper-paramètres (e.g., variance de la distribution normale de la proposition, ou variance du moment initial pour HMC).

Ces méthodes produisent des chaînes de valeurs de paramètres (échantillons). L’utilisation de tel ou tel algorithme MCMC pour échantillonner la distribution postérieure repose sur trois objectifs :

• Les valeurs de la chaîne doivent être représentatives de la distribution postérieure
  \(~\bullet~\) Ces valeurs ne doivent pas dépendre du point de départ
  \(~\bullet~\) Ces valeurs ne doivent pas être cantonnées à une région particulière de l’espace des paramètres
  

• La chaîne doit être suffisamment longue pour assurer la précision et la stabilité du résultat
  \(~\bullet~\) La tendance centrale et le HDI calculés à partir de la chaîne ne doivent pas changer si on relance la procédure
  

• La chaîne doit être générée de manière efficace (avec le moins d’itérations possible)

Évaluation des MCMC - Représentativité

Évaluation des MCMC - Représentativité

Cette affichage ne montre que les 500 premières itérations. Les trajectoires ne se superposent pas au début (zone orange). La densité est également affectée. En pratique on supprime ces premières itérations (burn-in ou warm-up).

Évaluation des MCMC - Représentativité

Vérification numérique des chaînes : Le shrink factor (aussi connu comme \(\hat{R}\) ou Rhat) est le rapport entre la variance inter-chaînes et intra-chaîne. Cette valeur devrait idéalement tendre vers 1 (on la considère comme acceptable jusqu’à 1.01).

Évaluation des MCMC - Stabilité et précicion

Plus la chaîne est longue et plus le résultat sera précis et stable. Si la chaîne s’attarde sur chaque position, et que le nombre d’itérations reste le même, alors on perd en précision. Il lui faudra plus d’itérations pour arriver au même niveau de précision. L’autocorrélation est la corrélation de la chaîne avec elle-même mais décalé de \(k\) itérations (lag).

Évaluation des MCMC - Stabilité et précicion

La fonction d’autocorrélation est représentée pour chaque chaîne (en haut à droite). Un autre résultat rend compte de la précision de l’échantillon : l’effective sample size, \(ESS = \frac{N}{1 + 2 \sum_k ACF(k)}\). Il représente la taille d’un échantillon non-autocorrélé extrait de la somme de toutes les chaînes. Pour une précision raisonnable du HDI, il est recommandé d’avoir un ESS supérieur à 1000.

Évaluation des MCMC - Stabilité et précicion

L’erreur standard d’un ensemble d’échantillons est donné : \(SE = SD / \sqrt{N}\). Plus \(N\) augmente plus l’erreur standard diminue. On peut généraliser cette idée aux chaînes de Markov : \(MCSE = SD / \sqrt{ESS}\). Pour une précision raisonnable de la tendance centrale, il faut que cette valeur soit faible.

Évaluation des MCMC - Implémentation via brms

library(tidyverse)
library(imsb)
library(brms)

d <- open_data(howell)
d2 <- d %>% filter(age >= 18)

priors <- c(
  set_prior("normal(150, 20)", class = "Intercept"),
  set_prior("normal(0, 10)", class = "b"),
  set_prior("exponential(0.01)", class = "sigma")
  )

mod1 <- brm(
  height ~ 1 + weight,
  prior = priors,
  family = gaussian(),
  data = d2, 
  chains = 4, # nombre de MCMCs
  iter = 2000, # nombre total d'itérations (par chaîne)
  warmup = 1000, # nombre d'itérations pour le warm-up
  thin = 1 # thinning (1 = no thinning)
  )

Évaluation des MCMC - Implémentation via brms

# combo can be hist, dens, dens_overlay, trace, trace_highlight...
# cf. https://mc-stan.org/bayesplot/reference/MCMC-overview.html

plot(
  x = mod1, combo = c("dens_overlay", "trace"),
  theme = theme_bw(base_size = 16, base_family = "Open Sans")
  )

Évaluation des MCMC - Implémentation via brms

library(bayesplot)
post <- posterior_samples(mod1, add_chain = TRUE)
post %>% mcmc_acf(pars = vars(b_Intercept:sigma), lags = 10)

Évaluation des MCMC - Implémentation via brms

summary(mod1)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + weight 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   113.90      1.92   110.20   117.68 1.00     3646     2735
weight        0.90      0.04     0.82     0.99 1.00     3751     2827

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     5.11      0.20     4.75     5.51 1.00     4056     2935

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Évaluation des MCMC - Implémentation via brms

Bulk-ESS fait référence à l’ESS calculé sur la distribution des échantillons normalisée par leur rang, et plus particulièrement autour de la position centrale de cette distribution (e.g., moyenne ou médiane). On recommande que le Bulk-ESS soit au moins 100 fois plus élevé que le nombre de chaînes (i.e., pour 4 chaînes, le Bulk-ESS devrait être d’au moins 400).

Évaluation des MCMC - Implémentation via brms

post %>% # rank plots
  mcmc_rank_overlay(pars = vars(b_Intercept:sigma) ) +
  labs(x = "Rang", y = "Fréquence") +
  coord_cartesian(ylim = c(25, NA) )

Résumé du cours

Nous avons introduit et discuté l’utilisation d’échantillonneurs (samplers) pour obtenir des échantillons issus de la distribution postérieure (non-normalisée). Ces échantillons peuvent ensuite être utilisés pour calculer différentes statistiques sur la distribution postérieure (e.g., moyenne, médiane, HDI).

Travaux pratiques

On s’intéresse à la performance économique des capitales rgdppc_2000 en fonction de deux paramètres : la rudesse du paysage (plus ou moins vallonné) rugged et son appartenance au continent africain cont_africa.

library(tidyverse)
library(imsb)

d <- open_data(rugged) %>% mutate(log_gdp = log(rgdppc_2000) )
df1 <- d[complete.cases(d$rgdppc_2000), ]
str(df1)

'data.frame':   170 obs. of  6 variables:
 $ isocode    : chr  "AGO" "ALB" "ARE" "ARG" ...
 $ country    : chr  "Angola" "Albania" "United Arab Emirates" "Argentina" ...
 $ rugged     : num  0.858 3.427 0.769 0.775 2.688 ...
 $ cont_africa: int  1 0 0 0 0 0 0 0 0 1 ...
 $ rgdppc_2000: num  1795 3703 20604 12174 2422 ...
 $ log_gdp    : num  7.49 8.22 9.93 9.41 7.79 ...

Travaux pratiques

Écrire le modèle qui prédit log_gdp en fonction de la rudesse du terrain, du continent, et de l’interaction de ces deux variables avec brms::brm(), en spécifiant vos propres priors. Examinez ensuite les estimations de ce modèle (interprétation, diagnostiques des MCMCs).

\[ \begin{aligned} \log(\text{gdp}_{i}) &\sim \mathrm{Normal}(\mu_{i}, \sigma_{i}) \\ \mu_{i} &= \alpha + \beta_{1} \cdot \text{rugged}_{i} + \beta_{2} \cdot \text{continent}_{i} + \beta_{3} \cdot (\text{rugged}_{i} \cdot \text{continent}_{i}) \\ \alpha &\sim \cdots \\ \beta_{1}, \beta_{2}, \beta_{3} &\sim \cdots \\ \end{aligned} \]

Proposition de réponse

priors2 <- c(
  set_prior("normal(0, 100)", class = "Intercept"),
  set_prior("normal(0, 10)", class = "b"),
  set_prior("exponential(0.01)", class = "sigma")
  )

mod2 <- brm(
  log_gdp ~ 1 + rugged * cont_africa,
  prior = priors2,
  family = gaussian(),
  data = df1,
  chains = 4, # nombre de MCMCs
  iter = 2000, # nombre total d'itérations (par chaîne)
  warmup = 1000 # nombre d'itérations pour le warm-up
  )

Proposition de réponse

summary(mod2)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: log_gdp ~ 1 + rugged * cont_africa 
   Data: df1 (Number of observations: 170) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              9.22      0.14     8.94     9.50 1.00     2690     2904
rugged                -0.20      0.08    -0.36    -0.05 1.00     2688     2830
cont_africa           -1.94      0.23    -2.39    -1.49 1.00     2463     2514
rugged:cont_africa     0.39      0.13     0.13     0.64 1.00     2430     2775

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.95      0.05     0.86     1.06 1.00     4156     2453

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Proposition de réponse

plot(
  x = mod2, combo = c("dens_overlay", "trace"), pars = "^b_",
  theme = theme_bw(base_size = 16, base_family = "Open Sans")
  )

Proposition de réponse

pairs(mod2)